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SUMMARY 


Two  different  methods  of  chromotomographic  reconstruction  applicable  to  the 
MWIR  were  investigated  by  the  Optical  Sciences  Center,  University  of  Arizona  for  The 
Air  Force  Research  Laboratory  (AFRI/SNHI)  at  Hanscom  Air  Force  Base,  Hanscom, 
MA.  The  first  method  involved  the  decomposition  of  the  datacube  into  principal 
component  space,  using  the  center  order  broadband  image  to  estimate  the  spatial 
principal  components,  then  using  this  estimate  to  deconvolve  an  estimate  of  the  spectral 
principal  components.  The  second  method  used  an  extension  of  the  traditional 
maximum-likelihood  approach  to  suppress  noise  in  reconstructions  of  simulated  low 
background  scenes. 


PREFACE 


This  final  report  describes  development  of  reconstruction  algorithms  for  MWIR 
chromotomography  by  the  University  of  Arizona’s  Optical  Sciences  Center  for  The  Air 
Force  Research  Laboratory  (AFRI/SNHI)  at  Hanscom  Air  Force  Base,  Hanscom,  MA. 
The  principal  investigator  was  E.  L.  Dereniak. 

The  authors  wish  to  thank  Tom  Hamilton  of  SMDTC,  Huntsville,  A1  for 
providing  the  512  x  512  InSb  snapshot  camera  used  to  collect  data  for  this  project. 
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I.  INTRODUCTION 


Chromotomographic  imaging  spectrometers  have  been  operating  successfully  in 
the  visible  region  of  the  electromagnetic  spectrum  for  several  years.  Recently,  however, 
there  has  been  heightened  interest  in  adapting  the  technology  to  wavelengths  beyond  one 
micron.  In  addition  to  the  group  at  Hanscom  AFB,  the  Optical  Detection  Lab  at  the 
Optical  Sciences  Center,  University  of  Arizona  has  developed  a  prototype  MWIR 
computed  tomography  imaging  spectrometer  (CTIS).  The  prototype  system  utilizes  a 
simple  binary  phase  disperser  fabricated  in  GaAs  and  a  512x512  InSb  snapshot  camera. 

Standard  reconstruction  techniques  such  as  Expectation-Maximization1  and 
MART2  have  been  used  with  some  success  to  reconstruct  object  cubes  from  the  raw 
snapshots  produced  by  this  instrument.  This  report  documents  investigations  of  two  new 
reconstruction  techniques.  The  first  method,  designated  here  as  “Principal  Component 
Method”,  shows  potential  for  increased  reconstruction  speed.  The  second  method, 
designated  here  as  “Extended  Maximum  Likelihood  Method”,  demonstrates  increased 
immunity  to  signal-independent  noise  in  the  raw  data. 
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II.  PRINCIPAL  COMPONENT  METHOD 

The  forward  action  of  tomographic  imaging  spectrometers  is  described  either  by 
theoretical  model3  or  by  experimental  measurement4.  The  goal  of  our  software 
development  is  to  invert  the  problem  to  form  an  estimate  of  the  object  cube  from  the  raw 
data  collected  by  the  spectrometer.  Traditional  techniques  such  as  pseudoinverse  and 
iterative  algorithms  have  been  applied  to  yield  acceptable  reconstructions.  More 
advanced  algorithms  are  required  to  improve  results  using  a  priori  information  about 
target  scenes  and  the  imaging  system. 

Experiments  performed  at  the  University  of  Arizona  apply  to  the  second- 
generation  MWIR  Computed  Tomography  Imaging  Spectrometer  (CHS).  The 
instrument  employs  a  simple,  binary  phase  Computer  Generated  Hologram  (CGH)  to 
produce  several  simultaneous  tomographic  projections  of  the  3D  object  cube  onto  a  large- 
format  Focal  plane  array.  The  result  is  a  non-scanning,  flash  imaging  spectrometer.  The 
current  MWIR  CTIS  operates  in  the  3.0-5.0  pm  band.  The  current  CGH  limits  the 
system  to  five  useful  projections  onto  a  512  x  512  focal  plane  array.  The  spectral 
accuracy  of  the  instrument  is  limited  to  about  21  bands  by  the  low  number  of  projections. 

Recent  work  has  investigated  the  decomposition  of  the  datacube  into  principal- 

component  space, 

f(x,  y,  x)= £  un  (*>  y)  •  v„  M  •  0) 

n=\ 

The  set  of  orthonormal  vectors  {vn}  represents  the  spectral  principal  components 
of  the  datacube  and  is  complete  in  the  spectral  space.  The  set  of  orthonormal  images 
{un}  are  the  principal  component  images  corresponding  to  {vn}.  The  values  of  w  mdicate 
the  contribution  of  the  n*  principal  component  to  the  composite  datacube.  R  represents 
the  number  of  non-zero  values  of  w ,  which  for  a  discretized  datacube  may  be  less  than  or 

equal  to  the  number  of  wavelength  samples. 

Any  order  in  the  CTIS  raw  image  can  be  expressed  analytically  using  the 

datacube  decomposition  from  (1).  Equation  (2)  expresses  the  raw  unage  orders  as  a 
convolution  of  the  principal  component  image  and  principal  component  spectrum. 
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(2) 


gj (x,y)=  £ W„  • : un  (x,  j>)*  *|y« (mj?  ■  Tj )•  tj (rrijf  •  ry  )g(r  •  ry )] 

n= 1 


The  order  index  is  j,  ry  is  the  unit  vector  expressing  the  diffraction  angle  of  the  order, 
and  m  is  a  scaling  factor  proportional  to  the  order’s  diffraction  angle.  v„(/nyr*ry)  is  a 
two-dimensional  representation  of  the  spectral  principal  component;  v„(myr -ry)=v„(^. 
Similarly,  fy[myrry]  is  the  2D  representation  of  the  CTIS  transmission,  responsivity, 

and  order  diffraction  efficiency  as  a  function  of  wavelength,  t(A).  This  equation  is 
derived  by  expanding  the  3D  tomography  expression  into  the  principal  components  basis 


set. 

The  potential  of  this  expression  for  the  raw  image  is  revealed  by  Fourier 
transforming  (2)  with  R=l,  p  =  ^  +  rjrj: 

Gj  fell)-  fell)'  K  (j>  ■ •  ij  )*  *?)  '  1 fy  )J  (3) 

This  expression  was  simplified  by  setting  the  ntj  term  to  1,  which  is  reasonable  for 
a  discretized  system  in  which  the  sampling  of  the  spectrum  corresponds  to  exactly  one 
pixel  of  dispersion  in  the  first  order.  The  Fourier  transform  of  an  order  is  the  product  of 
the  Fourier  transform  of  the  spatial  principal  component  and  the  2D  representation  of  the 
Fourier  transform  of  the  spectral  principal  component. 

To  form  an  estimate  of  the  datacube  using  this  technique,  it  is  necessary  to  impose 
a  constraint  that  the  datacube  is  dominated  by  the  first  principal  component  (wi  >  10w2). 
Several  non-overlapping  orders  must  also  exist  in  the  raw  image.  The  initial  estimate  for 
uo  is  the  center  order  from  the  raw  image.  Dividing  Gj  by  Uo  yields  an  estimate  of  the 

principle  component  spectrum.  Terms  ry  and  Tj  are  known  from  calibration  data,  as 

well  as  small  corrections  in  the  diffraction  angle  due  to  tilt  of  the  CGH. 

Figure  1  shows  the  extracted  order  images  from  an  MWIR  CTIS  raw  image  of  a 
warm  coffee  cup.  Figure  2  shows  the  Fourier  transforms  of  the  center  and  diffraction 
images  of  the  same  image.  Multiple  orders  contribute  more  estimates  of  vo,  however, 
some  orders  will  have  better  estimates  of  some  parts  of  the  spectrum  than  others  (for 
instance,  notice  the  structure  oriented  at  about  145°  in  the  center  order  of  Figure  2  is 
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represented  strongly  in  the  lower  left  and  upper  right  diffraction  image  spectra,  where  the 
diffraction  angle  is  nearly  orthogonal  to  the  structure).  Thus,  a  regularized  ratio  may  be 
used  to  reduce  noise  and  weight  the  contributions  to  the  spectra.  This  technique  assumes 

that  tj(j l)  is  constant  for  all  orders  used,  but  that  is  a  reasonable  assumption  for  a 

properly  designed  CGH. 


Figure  1.  Order  images  extracted  from  an  MWIR  CHS  raw  image  of  a  warm  coffee  cup.  The  image  at  left 
is  the  center  order.  The  arrangement  of  the  diffraction  images  corresponds  to  their  alignment  on  the  FPA. 
The  image  satisfies  the  requirement  of  a  dominant  first  principal  component  since  the  scene  has  nearly 
uniform  emissivity  and  temperature. 
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There  are  two  constraints  which  may  be  imposed  to  improve  the  spectrum  estimate.  The 
first  principal  component  spectrum  must  be  positive  and  it  must  be  zero  outside  of  the 
transmission  band  of  the  spectrometer.  Since  these  constraints  are  both  imposed  on  the 
real  spectrum,  an  iterative  retrieval  algorithm  may  be  used  to  recover  data  at  frequencies 
not  present  in  the  object.  The  results  of  an  example  recovery  using  only  the  positivity 
constraint  are  demonstrated  for  the  coffee  cup  image  in  Figure  3. 


Figure  3.  Results  of  an  iterative  retrieval  of  the  first 
principal  component  from  an  MWIR  image  of  a 
coffee  cup. 


Equation  3  may  be  used  to  estimate  the  first  principle  component  image,  «0,  from 
the  diffracted  orders  and  the  estimate  of  vo.  In  this  case,  dividing  Gj  by  Vo  yields  an 
estimate  of  the  principle  component  image.  The  estimation  of  the  full  (^,rj)  plane  is  far 
from  complete  for  the  second  generation  MWIR  CHS  since  the  raw  image  has  only  four 
orders  (and  only  two  independent  r  directions). 

As  with  the  principal  component  spectrum  estimate,  constraints  in  the  image 
space  allow  an  iterative  recovery  of  parts  of  the  missing  spatial  frequency  data.  The 
constraints  are  similar  to  those  of  the  spectrum  since  the  boundaries  of  the  field  stop  are 
known  and  the  first  component  must  be  non-negative  (subject  to  the  non-negativity 
constraint  of  the  image).  In  addition,  it  is  possible  to  impose  the  constraint  that  the 
component  images  must  be  zero  wherever  the  center  order  is  zero.  This  approach  would 
only  be  unreasonable  if  the  diffraction  order  efficiency  into  the  center  order  were  very 
low  for  any  wavelength  range,  which  is  not  hue  of  any  of  the  CGH  gratings  currently  in 
use. 
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The  results  of  a  spatial  component  estimate  are  shown  in  Figure  4.  Part  (a)  clearly 
shows  a  recovery  of  the  horizontal  and  vertical  structures  (corresponding  to  crossed  sine 
functions  present  in  the  FFT  of  a  square  cup  and  a  square  field  stop).  The  diagonal 
structures  are  an  artifact  of  the  limited  number  of  projections  and  are  clearly  evident  in 
the  corresponding  image  (Figure  4b). 


(a)  (b) 

Figure  4.  Advanced  initial  estimate  of  u0  for  the  coffee  cup  image  (Figure  1).  (a).  Center  image  and 
positivity  constraints  applied  to  U0(Z,,r\)  (logarithmic  scaling),  (b)  FFT  of  after  15  iterations  of  the 

constraint  algorithm. 

Reassembling  the  initial  estimates  of  «o(x,y)  and  v0(X)  into  a  CTIS  raw  image 
estimate  (g(x,y))  according  to  (2)  recovers  for  83%  of  the  image  data  (root-square  sum). 
About  5%  of  the  image  data  residual  can  be  attributed  to  noise.  Since  the  raw  image  does 
not  contain  the  complete  object  cube  information,  the  estimated  datacube  would  have  a 
lower  accuracy  compared  to  the  object  cube. 

Subtracting  the  estimate  from  the  raw  image  leaves  a  new  raw  image  without  the 
first  principal  component.  An  attempt  to  recover  the  second  component  is  possible,  but 
most  of  the  constraints  which  applied  to  the  first  principal  component  are  no  longer  valid. 
The  center  order  is  not  a  reasonable  estimate  of  the  second  principal  component  image 
and  the  positivity  constraint  for  either  the  spectrum  or  the  image  do  not  apply.  Only  the 
constraints  on  image  and  spectrum  boundaries  may  be  employed. 

Based  on  experiments  using  several  images,  the  second  estimate  typically 
recovers  an  additional  8%  of  the  raw  image  information.  Instead  of  estimating  the 
second  component,  the  first  component  may  be  used  to  form  a  datacube  which  is  used  as 
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an  initial  estimate  for  an  iterative  reconstruction.  With  more  orders  or  a  raw  image  with 
less  noise,  the  contribution  of  the  second  component  may  become  useful. 
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HI.  EXTENDED  MAXIMUM  LIKELIHOOD  METHOD 


Of  particular  interest  to  our  group  at  the  Optical  Detection  Lab  are  tomographic 
reconstruction  algorithms  suitable  for  MWIR  low  background  applications  due  to  our 
recent  involvement  with  BMDO  and  due  to  recent  interest  from  members  of  the 
astronomy  community.  The  problem  of  detecting  a  thermally  dim  object  against  a  dark 
background,  which  may  be  dominated  by  camera  system  noise  rather  than  photon  noise, 
is  one  that  requires  a  modification  of  our  conventional  approaches.  Heretofore  we  have 
assumed  that  either  signal  dependent  or  signal  independent  noise  alone  characterized  the 

entire  image. 

Recently,  a  step  toward  the  solution  of  the  mixed  noise  problem  was  undertaken 
by  Garcia  and  Dereniak5.  The  MERT  algorithm  achieved  some  success  in  reconstructing 
images  in  which  the  background  was  dominated  by  system  noise  and  the  “bright” 
portions  of  the  image  were  relatively  dim.  However,  system  noise  suppression  was  not 
entirely  satisfactory.  Furthermore,  a  literature  search  revealed  that  algorithms  proven  to 
deal  effectively  with  this  situation  were  generally  too  slow  for  our  purposes.  This  section 
details  the  derivation  and  initial  testing  of  a  new  algorithm  designed  to  remedy  these 
problems. 


Derivation 

The  derivation  begins  identically  to  MERT5,  but  with  one  critical  difference.  The 
signal  dependent  photon  noise  is  treated  as  a  normal  point  process  rather  than  a  normal 
continuous  process.  This  approach  is  an  inherently  more  accurate  description  of  photon 
statistics  arising  from  the  quantum  nature  of  light.  As  before, 


g  =  Hf  +  nl  +  n2  (5) 

The  elements  of  the  measurement  vector  g  represent  the  number  of  detected  photons  per 
pixel  during  one  integration  time  in  our  imaging  system;  the  elements  of  the  object  vector 
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f  represent  the  number  of  photons  emitted  during  one  integration  time  for  each  resolution 
element  of  the  object.  The  matrix  H  is  essentially  a  discrete  representation  of  the 
imaging  system’s  point  spread  function.  Both  vectors  are  convenient  ways  of  expressing 
a  two-dimensional  array  in  a  one-dimensional  format  to  ensure  that  H  can  be  expressed 
as  a  two-dimensional  matrix  rather  than  a  rank  four,  or  higher,  tensor.  The  matrix  element 
Hy  gives  the  contribution  of  object  element  fj  to  measurement  element  g\.  We  wish  to 
find  the  most  likely  noiseless  object  f  that  could  have  produced  the  measurement  g  in  the 
presence  of  photon  noise,  represented  by  nl,  and  post-detection  noise,  represented  by  n2. 

To  apply  the  Maximum  Likelihood  principle,  we  will  assume  that  the  signal- 
independent  system  noise,  n2,  is  a  Normal  continuous  distribution  with  zero  mean  and 
standard  deviation  as,  which  can  be  measured.  We  will  further  assume  that  each  Poisson 
distributed  (Hf)m  +  nlm  is  well  approximated  by  a  signal-dependent  Normal  point 
distribution  with  mean  (Hf)m  and  variance  (Hf)m.  This  holds  true  when  (Hf)m  is  greater 
than  or  equal  to  ten6.  With  these  assumptions,  gm  is  the  sum  of  two  statistically 
independent  normal  random  variables,  one  continuous  with  a  zero  mean  and  variance  gs 
and  the  other  discrete  with  mean  and  variance  equal  to  (Hf)m.  The  probability  density  of 
measuring  gm  photons  in  the  presence  of  photon  noise  at  the  m*  detector  given  some 
object  emission  distribution  f  is 


p,fejf)=£ 


-V  2*(Hf). 


reXPl 


[g.-(Hf)J2 

2(Hf)m 


(6) 


The  summation  of  weighted  delta  functions  in  this  expression  reflects  the  fact  that  the 
photon  noise  is  quantized.  The  probability  density  of  measuring  gm  photons  in  the 
presence  of  the  system  noise  floor  is 


Pli&m  I  °s)  = 


1 

*$2k<7s 


(7) 
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The  probability  density  in  the  presence  of  both  noise  sources  is  the  convolution  of  the 
two  densities. 


Prigm  \f><Ts)  =  Pligm  I  f)  *  Pl(gm  I  <* ,) 


(8) 


or  explicitly, 


liver. 


CO  00 

(V 

'  L-w-Fl 

►exp< 

Sm  1 

J2,exP< 

-00 1=0 

[  2(Hf).  J 

2^2  J 

(9) 


Carrying  out  the  integration  and  considering  all  detector  positions  as  statistically 
independent  gives  the  probability  of  seeing  some  image  g  for  an  object  f. 


Pt  (8  I  f  ) 


1 

litas 


M  wj  I 

nzH- 

m- 1  i»0 


[/  -  (Hf  )m  ]2 1 
2(Hf)n  j 


(10) 


Taking  the  natural  log  in  the  hope  of  producing  amore  tractable  function  of  the  object 
estimate  f  yields 


In 


\pT  (g  |  f ,  es  )]=  “  ln(2  nas )  +  £  In 

m=l 


(11) 


To  find  a  maximum  likelihood  estimate  of  the  object,  we  take  the  partial  derivative  with 
respect  to  a  test  estimate  of  the  number  of  photons  emitted  from  the  nm  element  and  insist 
that  it  equal  zero. 

4-ln[pr(g|f,<r,)]=0  (12) 

Sfn 
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Using  the  fact  that 


d(Hf), 

Vn 


dfn 


(13) 


we  arrive  at  an  equality  which  maximizes  the  likelihood. 


We  can  formulate  an  iterative  algorithm  for  the  best  current  estimate  of  each  element  of  f 
by  using  (10  )  above  in  a  correction  factor  which  is  applied  to  the  previous  estimate  of 
each  element  of  f. 


Implementation  of  this  algorithm,  as  written,  is  extremely  computationally  intensive, 
making  it  too  slow  for  our  purposes.  Even  though  the  summations  over  /  can  be  truncated 
to  plus  or  minus  5  as  about  gm  ,  if  as=  10  then  conceivably  as  many  as  M  x  200  terms 
must  be  computed  for  each  element  of  f  where  M  is  the  number  of  elements  of  g.  Surely 
there  must  be  a  faster  way. 

Fortunately  the  sums  over  /  appear  as  a  ratio,  raising  the  possibility  of 
approximating  the  numerator  and  denominator  by  the  area  under  the  continuous  curve 
defined  by  the  product  of  normal  functions  rather  than  explicitly  adding  up  the  value  of 
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this  product  at  many  discrete  points.  If  these  areas  can  be  expressed  in  closed  form,  then 
significant  computational  savings  will  result.  The  concept  is  illustrated  in  Figure  5. 


Figure  5.  A  ratio  of  two  envelopes  evaluated  at  many  discrete  points 
can  be  approximated  by  the  ratio  of  the  areas  under  the  envelopes 
which  hopefully  can  be  expressed  in  closed  form.  The  approximation 
should  be  good  when  the  spacing  of  the  discrete  evaluation  points  is 
gmall  enough  that  the  envelope  doesn’t  change  appreciably  in  between 
points. 

The  approximations  involved  can  be  expressed  by  the  following  relations. 
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and  the  reduction  to  closed-form  can  be  written 


The  integrals  have  been  carried  out  to  19  standard  deviations  past  the  mean,  (Hf)m ,  of 

the  first  exponential  term  in  the  integrands  to  approximate  an  infinite  upper  limit.  Little 
computational  penalty  is  incurred  by  such  mathematical  overkill  by  inspection  of  the 
explicit  forms  of  FI  and  F2.  The  expressions  themselves  do  not  become  more  complex 
as  the  upper  limit  increases  provided  it  remains  finite.  Only  coefficients  in  the 
expressions  change. 
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With  a  =  (Hf  *)m ,  b  =  gm ,  c  =  of ,  the  algorithm  can  be  written  concisely  as 


M  U 

V1  nmn 

'  F\ 

(Hfi)ni,gm,a7r 

^3 

K> 

(Hf  *)m,g.,o*l 

M 


(22) 


Since  the  error  function  can  be  easily  implemented  in  a  look-up  table,  FI  and  F2  can  be 
considered  general  closed-form  approximations  to  the  nasty  summations  over  /  in  (15). 

Discussion 

We  wish  to  note  in  passing  the  similarity  between  this  algorithm  and  algorithms 
developed  by  Snyder7  and  Llacer  and  Nunez8  to  deconvolve  the  initially  blurred  Hubble 
Space  Telescope  images.  Their  techniques  modeled  the  photon  noise  as  a  Poisson  point 
process  and  the  system  noise  as  a  Gaussian  continuous  process.  They  arrived  at  an 
expression  of  similar  computational  complexity  to  (15)  but  were  unable  to  find  general 
closed-form  approximations.  Consequently  these  methods  produced  superb  results,  but 
ran  with  agonizing  slowness  except  for  certain  special  cases  when  approximations  could 
be  made.  This  is  not  a  problem  if  an  overnight  run  to  reconstruct  a  single  image  is 
tolerable.  Initial  comparisons  of  runtimes  indicate  that  our  new  technique  is  roughly  30 
times  faster. 

Initial  tests  of  the  new  algorithm  against  its  predecessor,  MERT,  were  conducted 
to  gauge  relative  performance  on  a  simple  reconstruction  problem.  An  optical  system 
with  an  x-y  separable  Gaussian  PSF  and  system  noise  equal  to  signal  photon  noise  was 
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simulated  with  a  Monte  Carlo  program  that  computes  individual  photon  trajectories  to 
accurately  simulate  quantum  spatial  noise.  The  original  17  x  17  pixel  object  and  noisy 
blurred  17  x  17  pixel  image  are  shown  below. 


Figure  6.  Original  object  (left)  and  noisy  blurred  image  (right).  The 
Gaussian  PSF  was  5  pixels  in  diameter  to  3o.  The  system  noise  floor 
was  10  input  referred  photons  and  the  bright  pixel  signal  level  was  100 
photons  per  integration  period. 


The  maximum  value  (bright  pixels)  in  the  image  was  100  photons  per  integration  period; 
the  minimum  value  (dark  pixels)  was  zero  photons  per  integration  period.  The  system 
noise  floor  was  10  photons  per  integration  period.  Convergence  properties  were  studied 
using  a  uniform  object  of  50  photons  per  integration  period  as  a  starting  point  for  the 
images  and  plots  below. 
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Figure  7.  Second,  fourth,  six  and  eighth  iterations  of  the  new 
algorithm  (left)  and  MERT  (right).  The  new  algorithm  deals  more 
effectively  with  system  noise  but  produces  an  offset.  The  black  comers 
are  artifacts  from  a  slightly  incomplete  H  matrix. 
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* 

Figure  8.  Plots  of  the  error  vector  f  —  f  versus  the  element  index  (1- 
289)  for  the  second,  fourth,  sixth  and  eighth  iterations  of  the  new 
algorithm  (left)  and  MERT  (right).  Note  in  particular  the  offset 
pedestal  associated  with  the  new  algorithm.  The  starting  point  for  both 
algorithms  was  a  uniform  object  of  50  photons  per  integration  period 
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The  convergence  properties  of  the  new  algorithm  can  be  conveniently  studied  in  a  plot  of 
average  error  versus  iteration  number.  Here  we  define  average  error  as 

(23) 
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Figure  9.  Average  error  versus  iteration  number  for  the  new  algorithm 
(solid  line)  and  MERT  (dashed  line).  The  primary  source  of  error  in  the 
new  algorithm  is  the  offset  pedestal  which  is  very  close  to  as .  The 
primary  source  of  error  in  MERT  is  pixel  to  pixel  reconstruction 
inaccuracy. 

As  mentioned  in  the  figure  caption  above,  the  primary  source  of  error  in  the  new 
algorithm  is  an  offset  pedestal  which  is  suspiciously  close  to  cts  .  We  have  inspected  plots 
of  the  error  vector  (see  Fig.  4)  for  large  iteration  numbers  and  for  different  combinations 
of  objects  and  system  noise  levels  and  found  the  same  phenomenon.  This  implies  that  the 
new  algorithm  could  be  dramatically  improved  simply  by  subtracting  the  value  of  as 
from  the  final  image  and  setting  negative  pixel  values  equal  to  zero.  However,  it  would 


NEW  =  solid  MERT  =  dashed 


Iteration  Number 


be  preferable  to  track  down  the  precise  point  in  the  algorithm  where  the  offset  is 
introduced  and  then  fix  the  problem  in  a  more  rigorous  fashion. 

So  far  we  have  described  a  new  image  reconstruction  algorithm  which  deals 
effectively  with  both  photon  noise  and  camera  system  noise.  Good  reconstructions  have 
been  obtained  for  an  image  with  a  symmetrically  blurred  point  spread  function,  but  this 
image  is  not  typical  of  a  raw  focal  plane  image  from  a  chromotomographic  instrument. 
In  addition,  cumbersome  mathematical  expressions  complicate  the  reconstruction 
algorithm  itself.  A  final  concern  is  that  the  algorithm  produced  an  image-wide  offset  that 
might  have  resulted  from  approximating  a  ratio  of  infinite  sums  by  a  ratio  of  definite 
integrals  in  the  course  of  the  derivation.  The  following  section  addresses  these  three 
issues. 

Algorithm  Simplification 

The  expression  for  Fl(a,b,c)  and  F2(a,b,c),  equations  (20)  and  (21)  in  our  first 
progress  report  have  been  simplified  by  considering  the  case  of  infinite  upper  limits  of 
the  integrals  in  (18)  and  (19)  from  which  they  were  derived.  The  simplified  expressions 
for  FI  and  F2  become 
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with  a  =  (Hf  *)m ,  b  =  gm ,  c  =  a] . 

The  simplified  algorithm  was  tested  against  the  exact  algorithm  with  open-form  sums  in 
(15)  and  the  results  are  shown  below  in  Figure  10. 


Figure  10.  Average  Error  per  pixel,  as  defined  in  (23)  versus  number 
of  iterations.  The  offset  pedestal  has  not  been  subtracted  here.  The 
same  test  object  and  noise  conditions  as  in  the  Figure  6  was  used. 


The  summations  in  the  open  form  expression  were  truncated  to  a  region  of  plus  or  minus 
5os  about  gm  to  avoid  excessively  long  reconstruction  times.  As  can  be  seen  from  the 
figure  the  ratio-of-integrals  approximation  is  virtually  indistinguishable  from  the  exact 
ratio-of-sums  solution. 
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Image  Offset 

As  mentioned  previously,  another  concern  was  that  the  ratio-of-integrals 
approximation  introduced  an  image- wide  offset.  Figure  10  addresses  this  concern  as  well 
because  it  is  apparent  that  both  the  exact  and  simplified  algorithms  converge  to  the  same 
value,  which  is  just  slightly  larger  than  the  system  noise  floor  of  ten  photon  arrivals.  In 
Figure  11,  the  second  and  eighth  iterations  of  the  two  algorithms  demonstrate  that  the 
primary  source  of  error  is  indeed  offset. 


Figure  11.  Plots  of  the  error  vector,  f  —  f ,  versus  element  number  for 
the  second  (top)and  eighth  (bottom)  iterations  of  the  simplified 
algorithm  (left)  and  the  exact  algorithm  (right).  Pixel  to  pixel  error 
from  noise  has  been  reduced  to  an  offset  pedestal. 

Accordingly  we  have  concluded  that  the  offset  pedestal  is  inherent  in  the  algorithm  and  is 
not  a  product  of  the  simplifying  approximations.  Although  we  cannot  prove  formally  at 
this  point  that  the  algorithm  converges  to  the  specified  value  of  the  system  noise  standard 
deviation,  we  will  subtract  it  from  the  final  image  from  here  on.  The  justification  is  that 
the  pedestal,  which  is  a  systematic  error,  can  be  subtracted  because  it  is  common  to  all 
pixels  or  resolution  elements  in  the  reconstruction.  Tests  have  shown  that  the  pedestal  is 
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determined  by  the  value  of  system  noise  input  to  the  algorithm  regardless  of  the  actual 
system  noise  in  the  input  data. 


Realistic  Simulation 

To  determine,  whether  the  new  algorithm  is  actually  useful  for 
chomotomography,  an  idealized  CHS  system  was  simulated  and  the  algorithm’s 
performance  was  evaluated  against  one  of  our  standard  reconstruction  methods,  MART. 
The  simulation  mapped  a  17x17x10  (x,y,X)  object  cube  onto  a  lllxlll  focal  plane  in 
four  projections  plus  the  center  order.  For  simplicity,  the  ten  spectral  slices  of  the  object 
cube  produced  exactly  ten  pixels  of  lateral  dispersion  in  each  projection.  Spectrally 
uniform  diffraction  efficiency  was  assumed.  Figures  12,  13,  14  show  the  object  cube,  the 
noiseless  focal  plane  image  and  the  focal  plane  image  with  temporal  Poisson  noise  and 
Gaussian  system  noise  with  a  standard  deviation  of  12  counts,  respectively. 


Figure  12.  Object  cube  consisting  of  ten  spectral  slices.  Wavelength 
increases  left  to  right,  top  to  bottom.  The  exact  spectral  bands  are 
irrelevant  because  the  dispersion  is  defined  in  terms  of  displacement  in 
pixels  in  the  focal  plane  image. 


Figure  13.  Noiseless  focal  plane  image.  Each  spectral  image  of  Figure 
12  has  been  displaced  radially  by  one  pixel  with  respect  to  its  short 
wavelength  neighbor. 


Figure  14.  Noisy  focal  plane  image.  Poisson  noise  was  introduced  to 
the  image  of  Fig.  4  and  added  to  Gaussian  noise  with  a  standard 
deviation  of  12. 


Figure  15.  MART  reconstruction  after  eight  iterations.  Signal- 
independent  noise  has  propagated  through  the  reconstruction  process. 
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Figure  16.  New  algorithm  reconstruction  after  eight  iterations.  Signal- 
independent  noise  has  been  well  suppressed  in  all  but  the  two  end 
bands.  It  is  believed  that  constant  diffraction  efficiency  versus 
wavelength  in  the  CTIS  simulation  caused  very  low  SNR  in  regions  of 
the  projections  where  spectral  overlap  is  minimal...  namely  the  end 
bands. 

Average  error  per  pixel,  was  calculated  for  the  two  reconstructions,  both  of  which  used  a 
uniform  initial  estimate.  MART  achieved  a  value  of  5.4  arrivals/pixel  whereas  the  new 
algorithm  achieved  a  value  of  3.3  arrivals/pixel.  The  reconstruction  time  for  the  new 
algorithm  was,  surprisingly,  only  27  %  longer  than  MART.  It  should  be  noted  that  a 
similar  result  could  not  have  been  achieved  by  a  3  x  3  median  window  filter  applied  to 
either  the  focal  plane  image  or  the  final  reconstruction  since  single  pixel/resel  objects 
would  have  been  obliterated. 


IV.  CONCLUSION 


Two  reconstruction  algorithms  are  described  in  this  report.  The  first  “Principal 
Component  Method”,  has  already  proven  useful  in  practice  for  providing  better  intial 
estimates  of  the  object  cube  to  enhance  the  speed  and  accuracy  of  traditional  iterative 
reconstruction  methods  such  as  MART  and  EM.  The  second,  “Extended  Maximum 
Likelihood  Method”,  should  prove  useful  for  anticipated  low  background  observations  in 


the  MWIR. 
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